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Abstract 

We use a combination of numeric and analytic techniques to determine the 
ground state phase diagram of the Bose-Hubbard Hamiltonian with longer 
range repulsive interactions. At half filling one finds superfluidity and an 
insulating solid phase. Depending on the relative sizes of near-neighbor and 
next near-neighbor interactions, this solid either follows a checkerboard or 
a striped pattern. In neither case is there a coexistence with superfluidity. 
However upon doping "supersolid" phases appear with simultaneous diagonal 
and off-diagonal long range order. 
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The interplay between different types of order in correlated boson and fermion mod- 
els is a rich area of current research. One normally thinks of this interplay as taking the 
form of a competition: In superconducting materials, pairing can compete with structural 
phase transitions (A-15's), charge density wave (CDW) order (B&i_ x K x Bi0 3 ) , and antifer- 
romagnetic order (high T c oxides). This competition is reflected in theoretical models like 
the Hubbard Hamiltonian [l] where it is believed that superconductivity and antiferromag- 
netism occupy disjoint regions of the phase diagram, and in electron-phonon models like 
the Holstein Hamiltonian || where a Peierls-CDW state usurps pairing at large interaction 
stengths near half-filling. Similarly, in bosonic systems with on-site repulsion, superfluid- 
ity competes with an insulating Mott phase. ||f|] Nevertheless, it has been suggested that 
in 4 He the solid and superfluid order might coexist. |5|,|6|] While there are only somewhat 
indirect experimental indications for such a phase, "supersolids" have been extensively 
studied with approximate analytic approaches. p|-p~3|j There have been to date very few 
Monte Carlo simulations. |T4]JT^| 

In this Letter, we will explore the question of supersolid order within the context of the 
two-dimensional (2D) Bose-Hubbard Hamiltonian 

H = -tJ2( a l a j + a ) a i) ~ V U i + V n i + V l H U i n 3 + ^2 XI n i n k- i 1 ) 

(ij) i i (ij) ((ife» 

Here <2j is a boson destruction operator at site i, and rii = a|aj. The transfer integral t = 1 
sets the scale of the energy, and fi is the chemical potential. Vq,Vx, and Vz are on-site, 
near-neighbor, and next-near-neighbor boson-boson repulsions. In the hard-core limit, the 
Bose-Hubbard Hamiltonian maps onto the quantum spin-1/2 Hamiltonian 

H = ~tJ2(StSr + SfSr) + V x £ S?S* + V 2 £ S?S' k - H z £ 5f . (2) 
(a) (ij) {{ik}) 

The field H z = fi — 2V\ — 2Vz. Ordering of the density corresponds to finite wave vector Ising 

type order, whereas superfluidity is described by ferromagnetic ordering in the XY plane. 

The mean field phase diagram of Eq. 2 has been worked out by Matsuda and Tsuneto 
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T6| . By increasing the density from half filling the following phases can be expected: a 
Neel state, corresponding to a checkerboard Bose solid with an ordering vector fc* = (jr, 7r); 
a ferromagnetic phase, with the net moment M xy ^ and M z ^ 0, corresponding to a 
superfluid; and a fully polarized magnetic phase, where only M z ^ 0, corresponding to a 
Mott-insulator. As the solid and the superfluid phases possess different broken symmetries, 
it is expected that the transition between them is first order. In an alternative scenario it was 
proposed that instead there could be two distinct second order transitions, where the two 
order parameters vanish at separate points. []6|,|9|,|l0| In the regime between the two transitions 



both order parameters are non-zero, hence it has been termed a supersolid. |9|,|l2|,[n| The mean 
field analysis revealed that longer range forces (V 2 ) are needed to stabilize the supersolid, 
although recently it was claimed that this conclusion changes in the soft core case, and 



a supersolid phase exists with nearest neighbor interaction alone [|17]]. These mean field 
results need to be reevaluated in the light of recent studies on the related Heisenberg model 
with competing first and second neighbor couplings J\ and J 2 which reveal the possibility 
of additional phases: a collinear phase, with alternating lines of up and down spins, at large 
J%j J\ |T%|-pU|], and a disordered phase at intermediate values of J2/J1 [pi],P^|. 



In this paper we address these questions first by extended mean field and spin wave 
calculations, then by Quantum Simulations. 

Mean Field and Spin Wave Analysis. We extend the mean field analysis by (i) 
considering the formation of a collinear phase, with ordering vector = (ir, 0) or (0, 7r), and 
(ii) introducing an approximate soft core representation by allowing the spin length S to be a 
variational parameter and adding a term H constraint = V — l) 2 to the Hamiltonian. We 

expand the ground state energy around the superfluid phase, and consider the eigenvalues 
corresponding to small oscillations of the density and superfluid order parameter, a procedure 
equivalent to generating a Ginzburg-Landau energy. Off half filling we find three phases: a 
superfluid, and collinear- and Neel- supersolids. The phase boundary between the superfluid 
and the collinear supersolid is located at V 2 = 1 + <5 2 /[4 — <5 2 + 24/Vo], whereas the phase 
boundary to the Neel supersolid is at V 2 = V x - 1 - 25 2 /[4 - 5 2 + 32/Vq]. Here 6 = p-l/2,p 



being the density of the boson gas. The phase diagram for p = 0.53 and Vq = 7 is shown in 
Fig. 1. It displays Neel- and collinear supersolid, and superffuid phases. At half filling the 
supersolid phases vanish, and two insulating solids are direct neighbors to the superfluid. 
This result is independent of Vo, ie it is true both in the hard and soft core limits. 

The analyses of the spin wave fluctuations which exist in the literature P JT2| , |2"3]] are in 
disagreement. The spectrum has been found to be either linear (T^j or quadratic |||23| at 



the density (or field-) tuned solid-supersolid phase boundary. This dependence is crucial 
for numerical studies, as it determines the dynamical critical exponent z and thereby the 
appropriate finite size scaling of the lattice. 

To settle the issue, we determine the spin wave spectrum allowing for the formation 
of a superfluid, a Neel solid and a Neel supersolid. The calculations involve a Bogoliubov 



transformation of coupled density and phase modes. Details will appear elsewhere. [24]] In 
the Neel solid there are two excitation branches in a halved Brillouin zone. Both branches 
are gapped. In the superfluid there is a Goldstone mode of linear k dependence at small 
k, and a well developed minimum around (it, it). Taking the continuum limit identifies this 
with the roton part of the helium dispersion. Finally, in the supersolid phase one has a 
gapless linear mode, and a gapped one, again in a halved zone. 

To clarify the physics of the transitions we determine the dispersion at the phase bound- 
aries. At the supersolid-Neel solid transition the critical mode is at small k. For the generic 
case, H z ^ 0, (no particle-hole symmetry in the boson language), the linear mode softens 
into a quadratic one, yielding a dynamical critical exponent z = 2. The disappearence of 
this Goldstone mode signals the destruction of superfluidity. This value of z agrees with 



that of Chester || and Cheng ||23|| , but differs from that of Liu and Fisher [pL2|| , who obtain 
z = 1. 

At the generic superfluid-to-supersolid transition the critical mode is at k = k*: solidifi- 
cation is signalled by the roton mode touching zero. The rotons also stiffen from a quadratic 
to a linear minimum, hence z = 1. At half filling, i.e. with particle - hole symmetry, both 
above transitions have z = 1. In a recent Monte Carlo study |l5j the same value for z was 



used in choosing the lattice size to resolve both the superfluid-supersolid and supersolid- 
solid transitions, whereas we find z = 1 and z = 2 respectively. Finally at high fields, at 
the superfluid-to-Mott insulator transition the Goldstone mode softens out again, leading to 
z = 2, in agreement with earlier field theoretical predictions || and numerical simulations 

Quantum Simulations. In the remainder of this paper we will discuss the results of 
Quantum Monte Carlo (QMC) simulations of the Bose-Hubbard Hamiltonian performed on 
the Connection Machine CM5. We use the world line QMC method in which the 2D quantum 
partition function Z of the Bose-Hubbard Hamiltonian is rewritten as a path-integral over a 
classical occupation number field n(j, r) by discretizing the inverse temperature (3 = L t At. 
We work in the canonical ensemble, mostly near the special "half-filled" point p = N^/N = 
1/2. 

The presence of solid ordering will be demonstrated by measuring the equal time density- 
density correlations, and their Fourier transform, the structure factor, 

S&) = ^£e* r (n(?,T)n(7 +/>)). (3) 
ji 

Long range solid order in the thermodynamic limit is signaled by a linear growth of S(k*) 
with the number of lattice sites, N, at some ordering vector k*. In Eq. 1, V\ drives the 
formation of a checkerboard phase with fc* = (71*, 7r), where sites are alternately empty and 
occupied, an Ising type Neel antiferromagnet in the spin language. V% favors a striped phase 
where lines of occupied sites in either the x or y direction alternate with lines of empty sites. 
In this case the structure factor peaks at either k* = (0, 7r) or (vr,0). We will measure the 
superfluid density by looking at a topological property of the boson world lines, the winding 
number [g,|25|. 

One can determine the ground state phase diagram either by simulating lattices with 
large (3 or else by appropriately scaling L T oc L z with linear spatial lattice size L. The 
latter technique assumes foreknowledge of z which is later justified by appropriate scaling 
behavior, but has advantages in the precise determination of phase boundaries. We shall 
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use it for that purpose when required. However, as has been done extensively in simulations 
of both fermion |26[] and boson systems, we will primarily choose (3 large enough so that 



observables no longer change and we are assured of measuring ground state properties. 

Results at Half-Filling. Our first question is whether the supersolid phase can 
exist without a finite V 2 , as suggested by recent work |T5| , |r7|] . Fig. 2 shows the staggered 



structure factor S(tt,7t) and superfluid density p s as a function of near neighbor interaction 
strength V\ at Vo = 7 and V 2 = 0. We see a sharp transition in both quantities at V\ ~ 5. 
The raw data already strongly suggest that there is no supersolid phase intervening between 
superfluid and solid. We performed the appropriate finite size scaling analysis and found that 
the transition points differ by at most 0.5%, which we regard as statistically insignificant. 
We can now try to drive the supersolid by turning on the next near neighbor repulsion. It 
seems reasonable to do so near the transition found in Fig. 2. Therefore, in Fig. 3 we show a 
plot of S(tt, 7r), S(tt, 0), S(0, 7r), and p s at Vo = 7, V\ = 5.5 (i.e. just inside the solid phase), 
sweeping V 2 . We see that V 2 induces a nonzero value of p s (only a weak V 2 is needed since we 
started so close to the transition), while simultaneously destroying the checkerboard order. 
For V 2 > Vi/2 we enter the striped solid phase (S(tt, 0) ^ 0), and again p s vanishes. S(n, ir) 
also shows a small kink at this superfluid-striped solid transition. 

We have also explored a case of very weak core bosons where Vo = 3. At V 2 = a 
transition from SF to checkerboard solid occurs at V\ =6.7, with no supersolid in between. 
Nor does turning on V 2 help create one. Our conclusion is that, in agreement with our mean 



field theory, and contrary to what has recently been found, |15[ no supersolid phase exists 
at p = 1/2 in the Bose-Hubbard model. 

Supersolids in the Defect Phase. We turn next to the doped phase where 5 = 
p — 1/2 ^ 0. In Fig. 4 we show a plot of p s and S(ti, n) versus VI at V 2 = and 5 = 0.03. 
We see that a tail of nonzero p s persists beyond the point where the solid has formed. Fig. 4 
contains data for two lattices sizes at the same doping, demonstrating that the tail is not a 
finite size effect. p s drops considerably, but remains nonzero, as the supersolid is entered. 
Our picture is that Nb = N/2 bosons freeze into a solid, leaving only the remaining bosons 
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mobile. These then condense into a superfluid. This is borne out by studying the height of 
the tail as a function of 5. We find that the height is proportional to 5, indicating that only 
defect bosons make up the superfluid condensate within the solid. We have confirmed that 
S(7r, 71") scales linearly with lattice size, so that long range crystalline order is indeed present 
in the checkerboard supersolid phase. 

Fig. 4 exhibits the superfluid in coexistense with the checkerboard solid. One can also 
get superfluidity in a striped solid phase. This is illustrated in Fig. 5 where we show 
S(tt,0), S(0,tt) and p s at Vq — 7, V± — 5.5 and 5 = 0.06 as a function of Vi. A nonzero 
value of p s is now present in the checkerboard supersolid at small Vi- Increasing V2 melts 
this supersolid and there is a large increase in p s as all the bosons participate in the con- 
densate. At yet larger V 2 the striped supersolid emerges. Note that in Fig. 5 we separately 
plot the superfluid fraction in the x and y directions. In the checkerboard supersolid and 
in the pure superfluid we find p sx = p sy . However, in the striped supersolid this rotational 
symmetry is broken and p sx 7^ p sy . The symmetry is broken randomly in the different runs, 
and there is the expected correlation between p sx [p sy \, and which of S(tt, 0)[S'(0, tt)} is large. 
The one dimensional superfluid flows only down the appropriate channels left open by the 
striped solid phase. The boson wavefunction is localized in the orthogonal direction. In the 
checkerboard super-solid a similar symmetry breaking occurs, namely p sa 7^ p s b on the a/b 



sublattices. [12 



For our zero temperature quantum phase transition, the static periodic potential of solid 
bosons does not confine defects in the Neel supersolid. Their wave functions are still extended 
Bloch states, from which they can condense into a superfluid phase. A doped boson can 
move either through an intermediate state of double occupation of energy cost A = Vo, or 
through a move of two neighboring bosons of energy cost A = 2V\. The excess bosons have 



extended superfluid wavefunctions with a mass renormalized to m* = Jl + (A/16t) 2 /2£ 
from the bare m = l/2t. That bosons can move in the solid without double occupancy 
explains that the supersolid is present in the hard-core limit as well. Meanwhile, in the 
striped supersolid, the doped bosons move entirely freely along the lines of unoccupied sites. 



We find p s is larger by about a factor of two in the striped supersolid than in the Neel 
supersolid at the same doping. 

In conclusion, we studied the formation of supersolid phases in interacting boson systems 
by mean field, spin-wave and quantum Monte Carlo techniques. The results of mean field 
analysis qualitatively agree with the phase diagram, obtained by the quantum simulations. 
We have shown that a supersolid phase, instead of existing in some special and narrow 



window of parameter space, is a rather generic feature of the Bose-Hubbard model. f28 



Defects or interstitials introduced into the checkerboard and striped solid phases do not 
destroy the diagonal long range order, but rather bose condense into a superfluid. It is as 
yet unclear whether this scenario for supersolids is realized experimentally. There is one 
positive and numerous negative experiments studying the existence of a supersolid phase 
in bulk 4 He |J . Some thin film studies indicate the existence of superfluidity in incomplete 
layers on top of close-packed solid ones. || In this situation one might imagine that different 
layers are giving rise to the two types of order, rather than a single layer being both solid and 
superfluid. However, the connection between the theoretical and experimental situations is 
not entirely clear, since the precise relation of layer formation and multiple occupancy in 
the Bose-Hubbard model is still somewhat uncertain. J^J 

We acknowledge useful discussions with D. Arovas. This work was supported by Na- 
tional Science Foundation grant DMR 92-06023 and by Thinking Machines Corporation. 
A.P.K. gratefully acknowledges support through a habilitation scholarship of the Deutsche 
Forschungsgemeinschaft . 



S 



REFERENCES 

[1] The Hubbard Model, Rasetti and Tosatti, World Scientific (1992). 

[2] R.T. Scalettar, D.J. Scalapino, N.E. Bickers, Phys. Rev. B40, 197 (1989). 

[3] M.P.A. Fisher, P.B. Weichman, G. Grinstein, and D.S. Fisher, Phys. Rev. B40, 546 
(1989). 

[4] G.G. Batrouni, R.T. Scalettar, and G.T. Zimanyi, Phys. Rev. Lett. 65, 1765 (1990). 
[5] M.W. Meisel, Physica 178, 121 (1992), and references therein. 

[6] A.F. Andreev, "Quantum Crystals," in Progress in Low Temperature Physics, Vol. VIII, 
D.G. Brewer (ed), North Holland, Amsterdam, (1982). 

[7] G.A. Lengua and J.M. Goodkind, J. Low Temp. Phys. 79, 251 (1990). 

[8] A.F. Andreev and I.M. Lifshitz, Sov. Phys. JETP 29, 1107 (1969). 

[9] G. Chester, Phys. Rev. A2, 256 (1970). 

[10] A.J. Leggett, Phys. Rev. Lett. 25, 1543 (1970). 

[11] I.E. Dzyaloshinskii, PS. Kondratenko, and V.S. Levchenkov, Sov. Phys. JETP 35, 823 
(1972). 

[12] K.S. Liu and M.E. Fisher, J. Low Temp. Phys. 10, 655 (1973). 

[13] C. Bruder, R. Fazio, and G. Schon, Phys. Rev. B47, 342 (1993). 

[14] E.Y. Loh, Ph.D. thesis, University of California, Santa Barbara (1985). 

[15] A. van Otterlo and K.-H. Wagenblast, unpublished. 

[16] H. Matsuda and T. Tsuneto, Suppl. Prog. Theor. Phys. 46, 411 (1970). 

[17] E. Roddick and D.H. Stroud, unpublished. 

[18] R.R.P. Singh and R. Narayanan, Phys. Rev. Lett. 65, 1072 (1990). 

9 



[19] E. Dagotto and A. Moreo, Phys. Rev. Lett. 63, 2148 (1989). 

[20] P. Chandra, P. Coleman and A. Larkin, J. of Phys. - Cond. Mat. 2, 7933 (1990). 

[21] N. Read and S. Sachdev, Phys. Rev. Lett. 66, 1773 (1991). 

[22] F. Figuerido et al, Phys. Rev. B 41, 4619 (1990). 

[23] Y.C. Cheng, Phys. Rev. B 23, 157 (1981). 

[24] R.T. Scalettar, G.G. Batrouni, G.T. Zimanyi and A. Kampf, unpublished. 
[25] E.L. Pollock and D.M. Ceperley, Phys. Rev. B36, 8343 (1987). 

[26] J.E. Hirsch and S. Tang, Phys. Rev. Lett. 62, 591 (1989); S.R. White, D.J. Scalapino, 
R.L. Sugar, E.Y. Loh, J.E. Gubernatis, and R.T. Scalettar, Phys. Rev. B40, 506 (1989). 

[27] G.T. Zimanyi, R.T. Scalettar, P. Crowell, and G.G. Batrouni, unpublished. 

[28] Part of the distinction between a supersolid "window" and a broader supersolid phase 
results from the difference between the nature of the phase diagram in the N-V and 
fi-V planes. A sweep at constant chemical potential can cut into a solid phase, revealing 
a supersolid window, while a sweep at fixed density skirts the pure solid and remains in 
the supersolid phase. This occurs because of the existence of a gap in the solid phase, so 
that even a slightly doped system may have n substantially shifted from the undoped 
case. 



10 



Figure Captions 

Fig. 1: The mean field phase diagram in the soft-core case. 

Fig. 2: p s and S(ir, ir) are shown as a function of V\ at p — 1/2, Vq — 7, and Vi = 0. The 
lattice size is 8x8x16 and r — 1/4. There is a single transition between superfluid and solid. 

Fig. 3: S(ir, 7r), S(ir, 0), 5(0, n) and p s are shown as a function of V 2 at p = 1/2, V = 7, 
and V\ = 5.5. The lattice is 8x8x16 and r = 1/4. V2 first melts the checkerboard solid into 
a superfluid, and then causes a striped solid to form. p s is zero in the solid phases. 

Fig. 4: p s and £(71-, 7r) are shown as a function of V 1 at V^o = 7.0, V 2 = 0.0, r = 1/4, 
with the system doped to p = 0.53. A superfluid tail persists after the formation of the 
checkerboard solid. Open symbols are used for p s and closed symbols for S. Triangles are 
8x8x16 lattices and squares are 10x10x16 lattices. 

Fig. 5: p s , S(n,0) and S(0,n) are shown as a function of V2 at Vo = 7.0, Vi = 5.5, and 
the system doped to p s = 0.56. Here filled and open squares are for S(ir, 0) and S(0, it) 
respectively. Crosses and open triangles show p sx and p sy respectively. Both solid phases 
have nonzero p s . In the checkerboard supersolid the mobile dopants move in 2D (both 
p sx and p sy are nonzero), while here, in the striped phase the superfluid is confined to the 
appropriate channels between the lines of occupied sites. S(ir, n) is small throughout the 
parameter range and is not shown. 
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